/*
 * MAtoms.cpp
 *
 *  Created on: Aug 7, 2013
 *      Author: marchi
 */

#include "MAtoms.h"


template <typename T, typename U>
class create_map{
private:
    map<T, U> m_map;
public:
    create_map(const T& key, const U& val){m_map[key] = val;}

    create_map<T, U>& operator()(const T& key, const U& val) {
    	m_map[key] = val;
    	return *this;
    }

    operator map<T, U>(){return m_map;}
};

template <typename T>
class create_vector{
private:
    vector<T> m_vec;
public:
    create_vector(const T& data){m_vec.push_back(data);}

    create_vector<T>& operator()(const T& data) {
    	m_vec.push_back(data);
    	return *this;
    }

    operator vector<T>(){return m_vec;}
};

int MAtoms::calls=0;
vector<string> * MAtoms::ResList0=NULL;
vector<int> * MAtoms::ResIndx0=NULL;

map<string,double> MAtoms::MMass=create_map<string,double>
("H",1.008)("C",12.011)("N",14.0067)("O",15.9994)("S",32.06)("P",30.974)
("FE",55.847)("CA",40.080)("ZN",65.37)("NA",22.98977)("MG",24.305);
string MAtoms::MIons="CA NA SOD CAL ZN";

class PBCvect{
	static vector<Dvect> * nxyz;
	const int M;
public:
	PBCvect(): M(3){
		nxyz=new vector<Dvect>;
		for(int o=0;o<M;o++){
			double nx=static_cast<double>(o-1);
			for(int p=0;p<M;p++){
				double ny=static_cast<double>(p-1);
				for(int q=0;q<M;q++){
					double nz=static_cast<double>(q-1);
					nxyz->push_back(Dvect(nx,ny,nz));
				}
			}
		}
	}
	static vector<Dvect> & getVec(){return *nxyz;}
};
struct op_iComp {
  bool operator() (double i,double j) { return (i>j);}
};

struct op_sortResPDB{
	bool operator()(string x, string y){
		stringstream ssx(x.substr(21,5));
		stringstream ssy(y.substr(21,5));
		int nx, ny;
		ssx>>nx; ssy>>ny;
		return (nx < ny);};
};
struct op_sortAtnoPDB{
	bool operator()(string x, string y){
		stringstream ss1x(x.substr(6,5));
		stringstream ss1y(y.substr(6,5));
		int nx, ny;
		ss1x>>nx; ss1y>>ny;
		return (nx < ny);
	}
};


class Jacob{
	double ** iner;
	double ** ei;
	double * Im;
	static Dvect Im_s;
	void  eigsort(Dvect  & d, Matrix  * v=NULL){
		int k;
		int n=DIM;
		for (int i=0;i<n-1;i++) {
			double p=d[k=i];
			for (int j=i;j<n;j++)
				if (d[j] >= p) p=d[k=j];
			if (k != i) {
				d[k]=d[i];
				d[i]=p;
				if (v != NULL){
					for (int j=0;j<n;j++) {
						p=(*v)[j][i];
						(*v)[j][i]=(*v)[j][k];
						(*v)[j][k]=p;
					}
				}
			}
		}
	}
public:
	Jacob():  iner(NULL), ei(NULL),Im(NULL){};
	~Jacob(){
		for(int o=0;o<DIM;o++) delete [] iner[o];
		delete [] iner;
		for(int o=0;o<DIM;o++) delete [] ei[o];
		delete [] ei;
		delete [] Im;
	}
	void operator()(Matrix & iner0, Dvect & Im0, Matrix & ei0, int & nrot){
		if(!iner) {
			iner=new double * [DIM];
			for(int o=0;o<DIM;o++) iner[o]=new double [DIM];
		}
		if(!ei) {
			ei=new double * [DIM];
			for(int o=0;o<DIM;o++) ei[o]=new double [DIM];
		}
		if(!Im) Im=new double [DIM];
		for(int o=0;o<DIM;o++){
			for(int p=0;p<DIM;p++){
				iner[o][p]=iner0[o][p];
			}
		}
		jacobi(iner,DIM,Im,ei,&nrot);

		for(int o=0;o<DIM;o++){
			Im0[o]=Im[o];
			for(int p=0;p<DIM;p++){
				ei0[o][p]=ei[o][p];
			}
		}

		eigsort(Im0,&ei0);
	}

};
Dvect Jacob::Im_s=0.0;

vector<Dvect> * PBCvect::nxyz=NULL;

class CellComp{
	Dvect Xd;
	Matrix co;
public:
	CellComp(Dvect& x, Dvect & y, Matrix & c): Xd(x-y), co(c){};
	CellComp(Dvect& x, Matrix & c): Xd(x), co(c){};
	bool operator()(const Dvect & x, const Dvect & y){
		Dvect x1=Xd+x;
		Dvect x2=Xd+y;
		Dvect xc1=co*x1;
		Dvect xc2=co*x2;
		return xc1.Norm() < xc2.Norm();
	}
};
void MAtoms::Percolate(vector<vector<int> > & Mols){
	vector<string> resn=vector<string>(PDB.size());
	for(size_t o=0;o<resn.size();o++){resn[o]=PDB[o].resn;}
	Perco=new Percolation(SelRes,rd,resn);

}
Percolation & MAtoms::Percolate() {
	try{
		if(!Perco) throw string("Should initialize percolation. Abort.");
	} catch(const string & s){
		cout << s <<endl;
		exit(1);
	}
	vector<Dvect> v(x,x+nr);

	Perco->doContacts(v);
	Perco->gCluster();
	Perco->Accumulate();
	return *Perco;
}
MAtoms::MAtoms(): Rg_count(0), Perco(NULL){
	// TODO Auto-generated constructor stub
}
void MAtoms::initLists(TopolPDB & y,vector<string> & L){
	PDB=vector<PDBdata>(y.Size());

	for(size_t o=0;o<PDB.size();o++) PDB[o]=y[o];
	ResList0=new vector<string>(L);
	ResIndx0=new vector<int>;
	int ia=0;
	for(size_t o=0;o<PDB.size();o++){
		string sub1=PDB[o].atn;
		string sub2=PDB[o].resn;
		if(MIons.find(sub2) == string::npos) mass.push_back(MMass[sub1.substr(0,1)]);
		else {
			mass.push_back(MMass[sub1.substr(0,2)]);
		}
		if(!ResList0 || find(ResList0->begin(),ResList0->end(),sub2) != ResList0->end()) {
			ResIndx0->push_back(ia);
		}
		ia++;
	}
}

void MAtoms::Gyro(){
	vector<vector<int> > mCluster=Perco->getCluster();
	vector<vector<int> > mAtoms=Perco->getAtoms();
	Rg_i=vector<Gyration>(mCluster.size());
	Rg_stat=vector<Gyration>(mCluster.size());

	Gyration::setTime(time_c);
	for(size_t o=0;o<mCluster.size();o++){
		Rg_i[o]=0.0;
		Dvect cm=0.0,cmG=0.0;
		double unit_nmm2=1.0/(unit_nm*unit_nm);

		int ntot=0;
		double tmass=0.0;
		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			for(size_t q=0;q<mAtoms[n].size();q++){
				int i=mAtoms[n][q];
				for(int o1=0;o1<DIM;o1++) cm[o1]+=x[i][o1]*mass[i];
				for(int o1=0;o1<DIM;o1++) cmG[o1]+=x[i][o1];
				tmass+=mass[i];
				ntot++;
			}
		}
		cm/=tmass;
		cmG/=static_cast<int> (ntot);
		Dvect x0;
		Dvect Gm,Im,axis;
		Matrix Giner=0.0;
		Matrix ei,ei2=0.0;
		double MyRg=0;
		int oa=0;
		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			for(size_t q=0;q<mAtoms[n].size();q++){
				int nn=mAtoms[n][q];
				for(int o1=0;o1<DIM;o1++) x0[o1]=cmG[o1]-x[nn][o1];
				Giner+=x0%x0*unit_nmm2;
			}
		}
		Giner/=static_cast<double> (ntot);
		int nrot;
		Jacob()(Giner,Gm,ei,nrot);
		for(int o1=0;o1<DIM;o1++) MyRg+=Gm[o1];
		Matrix Inertia=0.0;
		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			for(size_t q=0;q<mAtoms[n].size();q++){
				int nn=mAtoms[n][q];
				for(int o1=0;o1<DIM;o1++) x0[o1]=cm[o1]-x[nn][o1];

				Inertia[XX][XX]+=(sqr(x0[YY])+sqr(x0[ZZ]))*mass[nn];
				Inertia[YY][YY]+=(sqr(x0[ZZ])+sqr(x0[XX]))*mass[nn];
				Inertia[ZZ][ZZ]+=(sqr(x0[XX])+sqr(x0[YY]))*mass[nn];
				Inertia[XX][YY]-=(x0[XX]*x0[YY])*mass[nn];
				Inertia[XX][ZZ]-=(x0[XX]*x0[ZZ])*mass[nn];
				Inertia[YY][ZZ]-=(x0[YY]*x0[ZZ])*mass[nn];
			}
		}
		Inertia[YY][XX]=Inertia[XX][YY];
		Inertia[ZZ][XX]=Inertia[XX][ZZ];
		Inertia[ZZ][YY]=Inertia[YY][ZZ];
		Jacob()(Inertia,Im,ei2,nrot);
		Im*=unit_nmm2;
		double fact=5.0/tmass/2.0;
		axis[XX]=fact*(Im[YY]+Im[ZZ]-Im[XX]);
		axis[YY]=fact*(Im[XX]+Im[ZZ]-Im[YY]);
		axis[ZZ]=fact*(Im[YY]+Im[XX]-Im[ZZ]);
		Rg_i[o](MyRg,Im,Gm,axis);
	}
	Rg_count++;
}
Dvect MAtoms::CM(){
	Dvect cm=0.0;
	double tmass=0.0;
	for(unsigned int n=0;n<ResIndx0->size();n++){
		int nn=(*ResIndx0)[n];
		for(int o=0;o<DIM;o++) cm[o]+=x[nn][o]*mass[nn];
		tmass+=mass[nn];
	}
	cm/=tmass;
	return cm;
}
Dvect MAtoms::CMCenter(){
	Dvect cm=0.0;
	double tmass=0.0;
	for(unsigned int n=0;n < MyRes.size();n++){
		for(unsigned int m=0;m < MyRes[n].size();m++){
			int nn=MyRes[n][m];
			for(int o=0;o<DIM;o++) cm[o]+=x[nn][o]*mass[nn];
			tmass+=mass[nn];
		}
	}
	cm/=tmass;
	return cm;
}
Matrix MAtoms::InertiaRotCenter(){
	Dvect x0,cm;
	cm=CMCenter();

	double unit_nmm2=1.0/(unit_nm*unit_nm);
	Matrix Inertia=0.0,ei,eiT;
	Dvect Im;
	int nrot;

	for(unsigned int n=0;n < MyRes.size();n++){
		for(unsigned int m=0;m < MyRes[n].size();m++){
			int nn=MyRes[n][m];
			for(int o=0;o<DIM;o++) x0[o]=x[nn][o]-cm[o];
			Inertia[XX][XX]+=(sqr(x0[YY])+sqr(x0[ZZ]))*mass[nn];
			Inertia[YY][YY]+=(sqr(x0[ZZ])+sqr(x0[XX]))*mass[nn];
			Inertia[ZZ][ZZ]+=(sqr(x0[XX])+sqr(x0[YY]))*mass[nn];
			Inertia[XX][YY]-=(x0[XX]*x0[YY])*mass[nn];
			Inertia[XX][ZZ]-=(x0[XX]*x0[ZZ])*mass[nn];
			Inertia[YY][ZZ]-=(x0[YY]*x0[ZZ])*mass[nn];
		}
	}
	Inertia[YY][XX]=Inertia[XX][YY];
	Inertia[ZZ][XX]=Inertia[XX][ZZ];
	Inertia[ZZ][YY]=Inertia[YY][ZZ];
	Jacob()(Inertia,Im,ei,nrot);
	eiT=ei.Transpose();
	return eiT;
}

Matrix MAtoms::InertiaRot(){
	Dvect cm=0.0,x0;
	double tmass=0.0;
	for(unsigned int n=0;n<ResIndx0->size();n++){
		int nn=(*ResIndx0)[n];
		for(int o=0;o<DIM;o++) cm[o]+=x[nn][o]*mass[nn];
		tmass+=mass[nn];
	}
	cm/=tmass;

	double unit_nmm2=1.0/(unit_nm*unit_nm);
	Matrix Inertia=0.0,ei,eiT;
	Dvect Im;
	int nrot;
	for(unsigned int n=0;n<ResIndx0->size();n++){
		int nn=(*ResIndx0)[n];
		for(int o=0;o<DIM;o++) x0[o]=x[nn][o]-cm[o];

		Inertia[XX][XX]+=(sqr(x0[YY])+sqr(x0[ZZ]))*mass[nn];
		Inertia[YY][YY]+=(sqr(x0[ZZ])+sqr(x0[XX]))*mass[nn];
		Inertia[ZZ][ZZ]+=(sqr(x0[XX])+sqr(x0[YY]))*mass[nn];
		Inertia[XX][YY]-=(x0[XX]*x0[YY])*mass[nn];
		Inertia[XX][ZZ]-=(x0[XX]*x0[ZZ])*mass[nn];
		Inertia[YY][ZZ]-=(x0[YY]*x0[ZZ])*mass[nn];
	}
	Inertia[YY][XX]=Inertia[XX][YY];
	Inertia[ZZ][XX]=Inertia[XX][ZZ];
	Inertia[ZZ][YY]=Inertia[YY][ZZ];
	Jacob()(Inertia,Im,ei,nrot);
	eiT=ei.Transpose();
	return eiT;
}
string MAtoms::gTotResList(){
	stringstream ss;
	for(vector<string>::iterator it=ResList0->begin();it != ResList0->end();it++)
		ss<< *it + " " ;
	return ss.str();
}
void MAtoms::InitSelection(vector<string> & y, TopolMic & MyTop){
	try{
		for(size_t o=0;o<y.size();o++){
			if(!MyTop.CheckResidue(y[o])) throw string(" Reference residue " + y[o] + " does not exist in the system. Abort");
		}
	} catch(const string & s){
		cout << s << endl;
		exit(1);
	}
	typedef map<string,vector<vector<int> > >  ResMap;
	typedef ResMap::iterator ResIter;

	ResMap & ResDef=MyTop.getDef();
	for(size_t o=0;o<y.size();o++){
		vector<vector<int> > Res=ResDef[y[o]];
		for(size_t p=0;p<Res.size();p++)
			SelRes.push_back(Res[p]);
	}
}
void MAtoms::Reconstruct(const string & y, TopolMic & MyTop){
	try{
		if(!MyTop.CheckResidue(y)) throw string(" Reference residue " + y + " does not exist in the system. Abort");
	} catch(const string & s){
		cout << s << endl;
		exit(1);
	}
	if(!calls++) {PBCvect dummy;}
	vector<Dvect> nxyz=PBCvect::getVec();
	typedef map<string,vector<vector<int> > >  ResMap;
	typedef ResMap::iterator ResIter;

	ResMap & ResDef=MyTop.getDef();
	vector<vector<int> > & Res=ResDef[y];
	MyRes=ResDef[y];
	Matrix co=Mt.getCO();
	Matrix oc=Mt.getOC();
	double v[DIM];
	v[XX]=co[XX][XX];v[YY]=co[YY][YY];v[ZZ]=co[ZZ][ZZ];
	double MinCO=*min_element(v,v+3)*0.5;  // Half the smallest axis is the cutoff distance for atom shifts
	Dvect Xref;

	// Apply boundary condition on every molecule as the default boundaries n .xtc files are atomic

	for(ResIter it=ResDef.begin(); it != ResDef.end();it++){
		vector<vector<int> > & Res0=it->second;
		for(unsigned int o=0;o<Res0.size();o++){
			Xref=xa[Res0[o][0]];
			for(unsigned int p=1;p<Res0[o].size();p++) {
				Dvect x1=Xref-xa[Res0[o][p]];
				Dvect xc1=co*x1;
				if(xc1.Norm() > MinCO) {
					Dvect tmp=*min_element(nxyz.begin(),nxyz.end(),CellComp(x1,co));
					for(int q=0; q < DIM;q++) xa[Res0[o][p]][q]=xa[Res0[o][p]][q]-tmp[q];
				}
			}
		}
	}
	// > Take the reference type molecules and eliminate boundary conditions --> begin

	// >>  Compute center of mass for each molecule


	vector<Dvect> xaa;
	for(unsigned int o=0; o< Res.size();o++){
		Dvect cm=0.0;
		for(unsigned int p=0; p< Res[o].size();p++){
			int pp=Res[o][p];
			cm+=xa[pp];
		}
		cm/=static_cast<double>(Res[o].size());
		xaa.push_back(cm);
	}


	// >> Obtain the center of mass for the cluster of the reference molecules

	Xref=xaa[0];
	for(unsigned n=1;n < xaa.size();n++){
		Dvect xb=xaa[n];
		Dvect tmp=*min_element(nxyz.begin(),nxyz.end(),CellComp(Xref,xb,co));
		xaa[n]=xaa[n]-tmp;
	}

	Dvect cma=0.0;
	for(unsigned n=0;n < xaa.size();n++){
		cma+=xaa[n];
	}
	cma/=static_cast<double> (xaa.size());
	// >>Shift all system by cma

	for(int n=0;n<nr;n++)
		for(int m=0;m<DIM;m++)
			xa[n][m]=xa[n][m]-cma[m];


	//> Apply molecular boundary conditions on the entire system
	for(ResIter it=ResDef.begin(); it != ResDef.end();it++){
		vector<vector<int> > & Res0=it->second;
		for(unsigned int o=0;o<Res0.size();o++){
			Dvect xcm=0.0;
			for(unsigned int p=0;p< Res0[o].size();p++) {
				xcm+=xa[Res0[o][p]];
			}

			xcm/=static_cast<double>(Res0[o].size());

			for(int m=0;m<DIM;m++) xcm[m]=rint(xcm[m]);
			for(unsigned int p=0;p<Res0[o].size();p++)
				for(int m=0;m<DIM;m++)
					xa[Res0[o][p]][m]=xa[Res0[o][p]][m]-xcm[m];

		}
	}

		//> Obtain new Cartesian coordinates from reduced xa's
	for(int i=0;i<nr;i++){
		for(int o=0;o<DIM;o++){
			xa[i][o]=xa[i][o]+0.5;
			x[i][o]=Mt.getCO()[o][XX]*xa[i][XX]+Mt.getCO()[o][YY]*xa[i][YY]+Mt.getCO()[o][ZZ]*xa[i][ZZ];
		}
	}
}
void MAtoms::Reconstruct(){
	if(!calls++) {PBCvect dummy;}
	vector<Dvect> nxyz=PBCvect::getVec();

	vector<vector<int> > mCluster=Perco->getCluster();
	vector<vector<int> > mAtoms=Perco->getAtoms();
	Matrix co=Mt.getCO();
	Matrix oc=Mt.getOC();
	double v[DIM];
	v[XX]=co[XX][XX];v[YY]=co[YY][YY];v[ZZ]=co[ZZ][ZZ];
	double MinCO=*min_element(v,v+3)*0.5;  // Half the smallest axis is the cutoff distance for atom shifts
	Dvect Xref;

	// Loop on the Clusters


	for(size_t o=0;o<mAtoms.size();o++){
		Xref=xa[mAtoms[o][0]];
		for(size_t p=0;p<mAtoms[o].size();p++){
			int n=mAtoms[o][p];
			Dvect x1=Xref-xa[n];
			Dvect xp=xa[n];
			Dvect xc1=co*x1;
			if(xc1.Norm() > MinCO) {
				Dvect tmp=*min_element(nxyz.begin(),nxyz.end(),CellComp(x1,co));
				for(int q=0; q < DIM;q++) xa[n][q]=xa[n][q]-tmp[q];
			}

		}
	}

	for(size_t o=0;o<mCluster.size();o++){
		Xref=0.0;

		vector<Dvect> xcm=vector<Dvect>(mCluster[o].size(),0.0);
		Dvect xcmC=0.0;

		// Compute molecule center of mass xcm

		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];

			// Apply boundary condition on every molecule as the default boundaries n .xtc files are atomic
			for(size_t i=0;i<mAtoms[n].size();i++){
				Dvect xx=xa[mAtoms[n][i]];
				xcm[p]+=xx;
			}
			xcm[p]/=static_cast<double>(mAtoms[n].size());
		}

		Xref=xcm[0];  // Reference to the first molecule of the cluster



		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			Dvect x1=Xref-xcm[p];
			Dvect xcmp=xcm[p];
			Dvect xc1=co*x1;

			if(xc1.Norm() > MinCO) {
				Dvect tmp=*min_element(nxyz.begin(),nxyz.end(),CellComp(x1,co));
				for(int q=0; q < DIM;q++) xcm[p][q]=xcm[p][q]-tmp[q];

				for(size_t i=0;i<mAtoms[n].size();i++){
					int ia=mAtoms[n][i];
					Dvect xaa=xa[ia];
					for(int q=0; q < DIM;q++) xa[ia][q]=xa[ia][q]-tmp[q];
				}
			}

		}
		// Compute the center of mass of the cluster

		for(size_t p=0;p<mCluster[o].size();p++){
			xcmC+=xcm[p];
		}
		xcmC/=static_cast<double>(mCluster[o].size());

		for(size_t p=0;p<mCluster[o].size();p++){
			int n=mCluster[o][p];
			for(size_t i=0;i<mAtoms[n].size();i++){
				int ia=mAtoms[n][i];
				xa[ia][XX]=xa[ia][XX]-rint(xcmC[XX]);
				xa[ia][YY]=xa[ia][YY]-rint(xcmC[YY]);
				xa[ia][ZZ]=xa[ia][ZZ]-rint(xcmC[ZZ]);
			}

		}

	}

	//> Obtain new Cartesian coordinates from reduced xa's
	for(size_t i1=0;i1<mCluster.size();i1++){
		for(size_t i2=0;i2<mCluster[i1].size();i2++){
			int n=mCluster[i1][i2];

			for(size_t i=0;i<mAtoms[n].size();i++){
				int ia=mAtoms[n][i];
				for(int o=0;o<DIM;o++){
					xa[ia][o]=xa[ia][o]+0.5;
					x[ia][o]=Mt.getCO()[o][XX]*xa[ia][XX]+Mt.getCO()[o][YY]*xa[ia][YY]+Mt.getCO()[o][ZZ]*xa[ia][ZZ];
				}
			}
		}
	}

}

MAtoms::~MAtoms() {
	// TODO Auto-generated destructor stub
}

ostream & operator<<(ostream & fout, MAtoms & y){
	y.cPrint(fout);
	return fout;
}

void MAtoms::cPrint(ostream & fout){
	vector<Dvect> xc=vector<Dvect>(this->nr);

	Matrix Rot0=InertiaRotCenter();
	Dvect cm=CMCenter();

	for(int n=0;n<nr;n++){
//		Dvect xb=x[n];
		xc[n]=x[n];
//		xb=xb-cm;
//		xc[n]=Rot0*xb;
	}


	stringstream ss;
	ss<<time_c;
	fout << string("REMARK    Generated by trjMicelles ") <<endl;
	fout << string("REMARK    SIMULATION TIME = "+ ss.str()) <<endl;

	vector<double> Par=Mt.getParas();
	ss.str(string());
	ss<< setw(9) << setprecision(3) << fixed << Par[0]*10.0;
	ss<< setw(9) << setprecision(3) << fixed << Par[1]*10.0 ;
	ss<< setw(9) << setprecision(3) << fixed << Par[2]*10.0 ;
	ss<< setw(7) << setprecision(2) << fixed << Par[3] ;
	ss<< setw(7) << setprecision(2) << fixed << Par[4] ;
	ss<< setw(7) << setprecision(2) << fixed << Par[5] ;
	fout<< string("CRYST1"+ ss.str() + " P 1           1") <<endl;

	map<int,vector<string> >MapRes;
	for(unsigned int n=0;n<ResIndx0->size();n++){
		string tmp;
		int o=(*ResIndx0)[n];
		ss.str(string());
		ss <<std::setw(5) << std::left << PDB[o].resn;
		string tmp2=PDB[o].first;
		tmp2.erase(26,1);
		tmp=tmp2.replace(21,5,ss.str());
		ss.str(string());
		ss << std::setw(5) << std::right << PDB[o].res+1;
		tmp.replace(21,5,ss.str());
		tmp.append(" 1.00");
		tmp.append(PDB[o].last);
		ss.str(string());
		ss << fixed << setw(8) << setprecision(3) << right << xc[o][XX]/unit_nm;;
		ss << fixed << setw(8) << setprecision(3) << right << xc[o][YY]/unit_nm;;
		ss << fixed << setw(8) << setprecision(3) << right << xc[o][ZZ]/unit_nm;;
		tmp.replace(30,24,ss.str());
		MapRes[PDB[o].res+1].push_back(tmp);
	}

	map<int,vector<string> >::iterator itt=MapRes.begin();
	for(;itt != MapRes.end(); ++itt){
		vector<string> & tmp=itt->second;
		std::sort(tmp.begin(),tmp.end(),op_sortAtnoPDB());
		for(size_t o=0; o < tmp.size(); o++){
			fout <<tmp[o] <<endl;
		}
	}
	fout << string("TER") <<endl;
	fout << string("ENDMFL") <<endl;
}
void MAtoms::fPrint(Fstream & fout){
	FstreamC * foutC=dynamic_cast<FstreamC *> (&fout);
	t_fileio *xd = foutC->getfin();
	gmx_bool bOK=bOK_c;
	int step=step_c;
	real prec=prec_c;
	real time=time_c;
	matrix box;

	vector<Dvect> x0=vector<Dvect>(ResIndx0->size());
	for(size_t i=0;i < x0.size();i++){
		int ia=ResIndx0->at(i);
		x0.push_back(x[ia]);
	}
	int    natos=static_cast<int> (x0.size());
	rvec * x=new rvec [natos];
	for(unsigned int m=0;m<x0.size();m++)
		for(int o=0;o<DIM;o++) x[m][o]=x0[m][o];

	for(int i=0;i<DIM;i++)
			for(int j=0;j<DIM;j++)
				box[i][j]=Mt.getCO()[i][j];
	try{
		if(!write_xtc(xd,natos,step,time,box,x,prec)) throw string("Cannot write next frame. Abort");
	}catch(const string & s){
		cout << s << endl;
		exit(1);
	}

}
Fstream & operator<<(Fstream & fout, MAtoms & y){
	y.fPrint(fout);
	return fout;
}


